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Abstract 

We present a new method to solve in a semianalytical way the Dokshitzer-Gribov-Lipatov- 
Altarelli-Parisi evolution equations at NLO order in the x-space. The method allows to construct 
an evolution operator expressed in form of a rapidly convergent series of matrices, depending 
only on the splitting functions. This operator, acting on a generic initial distribution, provides 
a very accurate solution in a short computer time (only a few hundredth of second). As an 
example, we apply the method, useful to solve a wide class of systems of integrodifferential 
equations, to the polarized parton distributions. 



PACS numbers: 12.38.-t, 13.60.-r, 02.60.Nm 



1 Introduction 



The scaling violation of nucleon structure functions is described in terms of Dokshitzer-Gribov- 
Lipatov-Altarelli-Parisi (DGLAP) evolution equations [y. The DGLAP integrodifferential equa- 
tions describe the Q 2 dependence of the structure functions, which are related, via the operator 
product expansion, to the parton distributions for which the DGLAP equations are usually written 
down. In this framework, the analysis of the experimental data, is performed fixing at some Qq 
the structure functions by assuming the parton distributions and computing the convolutions with 
the coefficient functions, which can be evaluated in perturbation theory. The comparison with 
experimental data, which are distributed at different values of Q 2 , goes through the solution of 
DGLAP equations for the parton distributions; thus a reliable and fast algorithm to solve these 
equations is welcome. 

In literature there are essentially three different approaches to solve the DGLAP equations. The 
first one is based on the Laguerre polynomials expansion This technique is quite accurate up 
to x- values not smaller than x ~ 10~ 3 ; on the contrary, below x the convergence of the expansion 
slows down M. Given that experimental data are already available down to about x, for the 
polarized case, and down to 10 -5 , for the unpolarized case, this method results no longer practical 
(for a more detailed discussion see section ^). 

An alternative approach takes advantage of the fact that the moments of the convolutions appearing 
in the equations factorize in such a way that the analytical solution, in the momentum space, can 
be written down |5|. However, also in the most favorable case in which the analytical expressions of 
the moments of the initial conditions are known, the numerical Mellin inversion is relatively CPU 
time consuming (see ||). 

Another strategy, the so called "brute force" method @, is fundamentally a finite-differences 
method of solution. It reaches a good precision in the small x-region |, H but requires a rather 
large amount of computer running time. 

In this paper we present a new semianalytical method, to solve DGLAP equations. It consists 
in constructing an evolution operator which, depending only on the splitting functions, can be 
worked out once for all. In this respect our strategy is similar to the one in [Q, ||. Our method 
to perform the convolutions instead, takes advantage of an x discretization (comparable to the one 
in ||, [|) which allows us to represent the evolution operator as a matrix. Thus the procedure to 
construct the solution reduces merely to a multiplication between the evolution matrix and an initial 
vector, and can be done in an extremely short computer time with the required accuracy. This 
is particularly appealing in the analysis of the experimental data on nucleon structure functions 
which requires a large number of parton evolutions. 

In the next section we discuss the (formal) analytical solution of the DGLAP equations; in the third 
one the algorithm to perform the x— integration is presented. The last two sections are devoted to 
analyze the numerical results relative to the evolution of polarized parton distributions, to study 
the yield of our method in comparison with others and to conclude. 
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2 The Evolution Operator 



In this section we review briefly the mathematical problem of DGLAP evolution equations and 
propose a method to solve them. In the following we limit ourselves to discuss the polarized parton 
distributions. The application of our method to the unpolarized case is straightforward. 

The DGLAP equation, up to Next-to-Leading-Order (NLO) corrections, for the Non-Singlet 
distribution is 



—Aq NS (x,t) = (APS(i) + a(t)AR NS (x)) <g> Aq NS (x,t) 
while for the Singlet and Gluon distributions we have: 

' APff(x) AP$(x) ' 



d_ / Aq s (x,t) 
8t { A~g(x,t) 

where 



APj?(s) AP<£\x) 



+ a(t) 



AR qq {x) AR qg {x) 
AR gq (x) AR gg (x) 



AR ij (x)^AP^\x)-^-AP^{x). 



2/?o 



In these equations the symbol (g) stands for 

f{x)®g{x)= 

Jx y \yj 



and 



f{x) = xf(x) 



Instead of Q 2 , we have used the variable t defined by 



t 



2 , 



a s (Q 2 



where a s is strong running coupling constant corrected at NLO: 



a s (Q 2 



4tt 



and so in Eqs. ([[])- (||) 



/3 /n(Q2/A 2 ) 
a(t) 



1 



fa ln(ln(Q 2 /A 2 )) 



<* s {Ql) 



2tt 



Exp 



{44 ^ 



(i) 



Aq s (x,t) 
Ag(x,t) 



(2) 
(3) 

(4) 
(5) 

(6) 

(7) 
(8) 



The explicit expressions for (3q, (5\ as well as for the Splitting Functions APij(x) can be found in 
0- 



The equations @j and @) can be written in the following general form: 





Of 



f(t) = fi(t) f (t) 



where f(i) indicates the "vector of components f(x,t)" and Q(t) a linear operator acting as: 



im ©f(t)L 



dy ui(x,y,t) f(y,t) . 



(9) 



(10) 
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Note that, in the Singlet-Gluon case, f(t) becomes a doublet of vectors and fi(t) a 2 x 2 matrix of 
operators. 

Due to the logarithmic dependence of t on Q 2 , the range of values of physical interest for t — to 
(to is the starting values of t, where the parton distributions are assumed known) is small enough 
to expect that the Taylor's series of the solution f(i) converges rapidly. On the other hand, by 
deriving repeatedly the Eq. lM) we can write: 



Qk 



f(t) 



f (t 



fill 



t=t 



where the operators M( fe ) can be obtained recursively: 

M<°> = I 

m (1) = n 



M< 2 > 
M< 3 > 

M (fc) 



(i) 



n 2) + 2 n 1} m (1) + n m (2) 



fc-i 

i=0 



(12) 



(k) 

The cj indicates the i— th term of the fe— th row of Tartaglia triangle and 



t=t 



n = n(to) , 

Then the solution can be written as: 

f < ' > = 1 E ) f (*o) = T(t - to) f (to) 



" (t-to) fc 

Vfc=0 



(13) 



(14) 



with T(t — to) the Evolution Operator. 

As we will point out in section [| the series in Eq. ( [l4| ) converge quickly enough to obtain a 
very good approximation retaining only a first few terms. It is worth to note that if the operator 
f2(t) can be written as h(t) Q,' (with h(t) a numerical function) it is easy to show that the series in 
Eq. (|i~4|) reduces to: 



f (t) = Exp I 



t 



h(r)dT 



& f (to) 



(15) 



This is the case of DGLAP equation at Leading Order (LO) approximation. Nevertheless, in 
Eqs. (|[)-(§), where NLO corrections are included, we have Q(t) = Q\ + a(t)fl2 , with Sl± and Q,2 
non-commuting operators. As a consequence the series in Eq. (|i~4|) cannot be summed and it is not 
possible write the solution in a closed form. 

In the next section we show how to represent the operators M( fc ) in a form suitable for numerical 
calculations. 
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3 The x- Integration 



The integrals in Eq. (||) are evaluated with a method that generalizes the one proposed in Ref. 
|J. The method consists to treat exactly the "bad" behaviour of the kernel u(x,y,t) in Eq. (10) 
and approximate the "smooth" function f(y,t). In particular, we construct a M + 1 points grid 
(xo > 0,xi, ...,xm-i,xm = 1) hi the interval ]0, 1] and approximate f(x) (here we omit for brevity 
the t dependence) in each subinterval as linear combination of a suitable set of function gi(x) (we 
have used the power series x 1 ^ 1 ): 

m 

f(x) fs a i 9i(%) Vx e [x k ,x k+1 ] . (16) 

1=1 

Here a^ are coefficients which depend linearly on the values f(x k ) of f(x) computed in the x k 
points; they can be obtained solving \/k the linear system of m equations, 

m 

f{Xk+v-m/2) = Q l 9l{Xk+r-m/2) r 6 {1, -, m } (17) 

1=1 

giving 

m 

«S fc) = E4 fe) /(x fc+ ,_ m/2 ). (is) 

r=l 

In our numerical calculations we choose m = 4 which corresponds to approximate f{x) in each 
interval [xi-, Xfe+i] by the cubic which fits the four point f(xi), with i = k — 1, k, k + 1, k + 2. In the 
subintervals near the end-points (xq, and 1) we have reduced the order m. 

The general structure of the Polarized Splitting Functions which appear in Eqs. ([l])-(|2]) isQ: 

AP(x) = .^ X \ + B{x) + 6(1 - x)C , (19) 
(l-x) + 

where 



lo (1 - z) + Jo 1 - z 

and therefore the "i component" of the convolution is: 

ap M «/M = „ ( /' * *«/»)/(■>) - *D/M + /' j| B (s) m ) + 

V./^ y y - Xi J Xi y l \y J J 

{C + A{l)ln{l-Xi))f{xi). (21) 

Substituting Eq. (Jl|) in Eq. @ we obtain Mi £ {0, M - 1} ( E^Lm = is understood) 

m Af— 1 m 

i=i fc=t+i i=i 

+ (c + ^(l)/n(l-x i )-^(l)a 4 )/(x l ) ; (22) 



"The same structure, however, holds for unpolarized and transversely polarized splitting functions. 
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then Eq. (Qq) implies: 



M 



AP(xi) <g> f(xi) = u ik f(x k ) 

k=0 



(23) 



where u is the matrix of the coefficients of f(x k ). The matrices /?, 7, p and cr are given in the 
Appendix. 



Due to the ln(l — xi) the Eq. (21) contains a logarithmic divergence for i = M. This is not a 
problem if the solution f(x,t) decreases to in x = 1 faster than the logarithm: in this case it's 
equivalent to put uJuk = Vfc in Eq. (|2^), but in general a numerical regularizator should be used 
to control the lim^^i ln(l — x)f(x,t). 



Therefore the Eqs. ([l|)-@ became [} 



0_ 

dt 



M 

Aq NS (xi,t) = ^ ( u ik NS + a (*Mfc NS ) AqNs{xk,t) 



(24) 



_5 / Aq s (xi,t) 
dt { Ag( Xi ,t) 



m r 

E 

fc=0 L 



k=0 



(0) 9 q , ,(0) qg 



ik 



(0) gq , ,(0) 99 



+ a(t) 



(1) 99 , ,W 99 



ifc 



(1) 99 , ,(1) 99 
ifc ^ik 



Aqs{x k ,t) 
Ag(x k ,t) 



(25) 

We solve these equations by means of the method shown in section |2|: the operator Q(t) and then 
the M( fc ) became now numerical matrices, and the symbol stands for the usual rows by columns 
product. We would stress the fact that the matrices M( fc ) depend only on the points Xi and so they 
can be numerically evaluated once for all. 



4 Numerical Analysis 



The convergence of our algorithm is controlled by two parameters: the order n of the truncated 
series 

" /< , J. 

-M^ , (26) 



T^(t-t ) = ± {t - t0)k - 



k=0 



which define the evolution operator, and the number M of the points of x— integration. 

To test the accuracy of our method we evolve the Gehrmann and Stirling polarized singlet-gluon 



xAu v + xAd v + 6xAS 
.918* 1.365 *x a512 (l-x) 3 - 96 



initial distributions (cf pl[| ): 
Aqs(x,t ) = 



.339 * 3.849 *x- 78 (l-x) 4 ' 96 



Ag(x,t ) 



(1 + 11.65a; - 4.6^) 
(1 + 7.81x - 3.48v^) 
6 * .06 * 18.521 * x' 724 (l - x) 14 - 4 (l + 4.63x - 4.96^) 
1.71 * 3.099 * x' 724 (l - x) 5 ' 71 



(27) 
(28) 



^Note that ur 1 ' matrices correspond to the convolutions of the parton distributions with AR (cf Eq. (0)). 
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from Ql = 4 GeV 2 (t = 0) to Q 2 = 200 GeV 2 (t = 0.136) and Q 2 = 50000 GeV 2 (t = 0.245). We 
choose to work, as in the paper Q, in the fixed flavour scheme, rif = 3, with Aqq D = 231 MeV, 
and without taking into account, in the Q 2 evolution of a s , quark thresholds. The range ]0, 1] has 
been divided in M steps by M + 1 points: xq,x\, xm distributed in such a way that the function 
ln(x) + 2x varies by the same amount at any step; this function is slightly different from the pure 
logarithmic distribution commonly used in literature ||, 0], but allow, in our case, a more uniform 
distribution of the numerical errors. The end points are fixed to be xo = 1 x 10~ 8 and xm = lj 
however, for a better reading, in the Figures [I]— ^ the x— axis ranges from 10~ 4 to 1. 

First, we fix M = 100. In Figs, [l]— g are reported the evolved singlet and gluon distributions, 
respectively, obtained with n = 3, 6 and 12 for Q 2 = 200 GeV 2 and Q 2 = 50000 GeV 2 . It is 
worth to note the very fast convergence of the series to the solution, as already observed above. 
As a matter of fact, the maximum difference between the solutions relative to n = 6 and n = 12 is 
1.4 x 10~ 5 (7.6 x 10" 4 ) for the singlet, and 9.3 x 10 5 (5.3 x 10 3 ) for the gluon distribution, in 
correspondence of Q 2 = 200 GeV 2 (Q 2 = 50000 GeV 2 ). 

Next we fix n = 12 and Q 2 = 200 GeV 2 . In Figs. |3|— ^ are plotted the approximated evolved 
distributions with M = 25, 50, 100: the maximum difference on the common points between 
M = 50 and M = 100 is 4.6 x 10 -4 for the singlet and 7.9 x 10 -4 for the gluons. By comparing 
the results in Fig. ^— § with the corresponding Figs. 1-4 in Ref. we observe, besides a good 
numerical agreement of the results, a faster convergence as the number M of integration points 
increases, as a consequence of our more accurate x-integration procedure with respect to the so 
called "brute force" methods. In fact it should be observed that reducing from m = 4 (cubic) to 
m = 2 (linear) the order of approximation of f(x) in Eq. ([!(]), the accuracy £(x,t) (defined in the 
sequel) becomes about 1 and 3 order of magnitude bigger, respectively for singlet and gluon. 

To discuss the degree of accuracy of our method, let us introduce a global accuracy £(x,t) 
defined as the difference between left and right-hand side of the Eq. (g). The comparison between 
the range of values of £ with the one of both sides of Eq. (||) represents a very good estimate of 
the degree of accuracy of the solution. In Figs. |5| and [6] are plotted, for n = 12, M = 100 and 
Q 2 = 200 GeV 2 both sides of the Eq. @ and the corresponding (rescaled) accuracy £(x,t) []. It 
appears evident that an excellent approximation of the solution is obtained. 

Another advantage of our method, once fixed the accuracy of the solution, appears to be the 
running time to get each evolution. In fact, the simple analytical structure of the evolution matrix 
T makes the solution procedure considerably fast. As a matter of fact, once given the Splitting 
Functions and constructed the corresponding matrices M( fc ) (we have used Mathematica []l2| to 
do this), a single evolution, i.e. the multiplication of the T evolution matrix by the initial vector, 
require, for n = 12 and M = 100, about 6 x 10~ 2 sec on an AlphaServer 1000 using a Fortran 
Code. 

''Note that the integration in the right-hand side has been performed numerically after an x-interpolation of the 
discrete values obtained with the evolution operator, while the left-hand side is worked out by direct derivation of 
Eq. (0). 
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Particularly interesting is the comparison between our method and the one presented in |2], |3| , 
where an evolution operator is also introduced. Firstly we observe that the latter method is based 
on a polynomial expansion of the splitting and distribution functions. The expansion is equivalent 
to an expansion in power of x. As a consequence it is affected by problems of convergence for 
x — > 0, due to the branch point in zero of the involved functions. This is the source of the difficult 
encountered in the small- a; region, which are not present in our approach in which an optimized 
Newton-Cotes-like quadrature formula is employed. 

Second, also in the x-region of convergence, the Laguerre polynomial expansion needs, for each 
evolution process, the computation of the moments of the initial parton distributions with respect 
to the polynomials: this procedure requires a remarkable amount of CPU-time with respect to our 
approach in which only the evaluation of the initial parton distribution in the M grid points is 
needed. 

5 Conclusions 

We have proposed a new algorithm to solve the DGLAP evolution equations, in the x-space, which 
appears suitable for a rather large class of coupled integrodifferential equations. 

The method produces a solution which is analytical in the Q 2 -evolution parameter and approx- 
imate, but rapidly convergent, in the x— space. It allows to construct, once for all, an evolution 
operator in matrix form. It depends only on the splitting functions appearing in the equations and 
can be rapidly applied to whatever initial distribution to furnish the evolved one, requiring for each 
evolution only a few hundredth of second. 

It is worth to note the reliability of our x— integration algorithm, which gets excellent approx- 
imations on the whole x— range (we use, for all the calculations, 10 -8 < x < 1), also with few 
integration points, resulting in an evolution matrix of particularly small dimensions. 

In conclusion, our method, whose numerical implementation is straightforward, appears to 
be very fast, very accurate and extremely stable with respect to the increasing of convergence 
parameters (i.e. n, the order of the truncated series which gives the Evolution Operator, and M, 
the number of integration points). For these reasons it represents a powerful tool to analyze the 
experimental data on nucleon structure functions. 
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Appendix 



Here we give the definitions j3, 7, p and a matrices appearing in Eq. ( p2| ) 

dy A{xj/y)gi(y) - A(l)gi(xj) 

y y - Xi 

d ^B{x, l /y)g l {y) (A.2) 

dy A{xj/y)gi{y) 
y y -Xi 
dy 











— ^7 


/Xi 


w 


— 






— 37? 








f - 



(A.4) 

where A, B and C are the coefficients appearing in the Eq. (|l9|). Note that, choosing ^(x) = x l_1 
it is possible to analytically evaluate all the previous integrals in term of logarithms and Li n (z) = 
X^Li z J /j n ; the final expressions are too long and are not reported here. 
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Figure 1: The initial Singlet distribution (Q 2 = 4 GeV 2 , solid line) and the evolved ones for n = 3 
(dashed lines), n = 6 (dotted lines) and n = 12 (solid lines) corresponding at Q 2 = 200 GeV 2 and 
Q 2 = 50000 GeV 2 . We use M = 100. 
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Figure 3: The initial Singlet distribution (Q 2 = 4 GeV 2 , solid line) and the evolved one at 
Q 2 = 200 GeV 2 with M = 100 (solid line), M = 50 (dotted line) and M = 25 (dashed line) with 
n= 12. 
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Figure 5: For n = 12, M = 100 and Q 2 = 200 GeV 2 both sides of the Eq. (@) are plotted (solid line 
and dotted line), in correspondence of the Singlet distribution. Dashed line represents £{x) x 10 3 
(see text). 
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